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We use the Bogohubov theory of Bose-Einstein condensation to study the properties of dipolar 
particles (atoms or molecules) confined in a uniform two-dimensional geometry at zero tempera- 
ture. We find equilibrium solutions to the dipolar Gross-Pitaevskii equation and the Bogoliubov-de 
Gennes equations. Using these solutions we study the effects of quantum fluctuations in the system, 
particularly focussing on the instability point, where the roton feature in the excitation spectrum 
touches zero. Specifically, we look at the behaviour of the noncondensate density, the phase fluctua- 
tions, and the density fluctuations in the system. Near the instability, the density-density correlation 
function shows a particularly intriguing oscillatory behaviour. Higher order correlation functions 
display a distinct hexagonal lattice pattern formation, demonstrating how an observation of broken 
symmetry can emerge from a translationally symmetric quantum state. 



I. INTRODUCTION 

The field of ultra-cold atoms is often heralded as an 
ideal system to solve outstanding problems in many-body 
physics, owing to the innate tunability of certain cru- 
cial parameters within the Hamiltonian. The extremely 
low densities and temperatures associated with such sys- 
tems typically ensure that interactions between atoms 
are completely dominated by 5- wave collisions. Tech- 
niques for varying the strength of such interactions, via 
Feshbach resonances, are long-standing [H, Exciting 
recent developments involve the introduction of an en- 
tirely new class of interaction, the dipole-dipole inter- 
action, by cooling atoms with a large magnetic dipole 
moment, such as chromium or dysprosium, to tempera- 
tures below quantum degeneracy [3-8]. Even more ambi- 
tious, are the attempts to cool heteronuclear molecules, 
which typically involve considerably larger dipole mo- 
ments [3, [loj. Such systems typically have permanent 
electric dipole moments. The energy associated with the 
dipole interactions of heteronuclear molecules can be ex- 
tremely large, greatly exceeding the thermal kinetic en- 
ergy, or any other type of interactions in the system. 
These exciting experiments suggest that the observation 
of strongly correlated regimes in degenerate dipolar gases 
are close at hand. Another recent experiment, also aimed 
at exploring many-body physics beyond s-wave interac- 
tions, demonstrated the feasibility of introducing angular 
dependency into the interactions by dressing the atomic 
states with light effectively generating screened in- 
teractions analagous to electron fluids in metals. These 
experiments and others of a similar vein, motivate the in- 
vestigation into the novel physics which potentially could 
be observed in systems with long-range and angularly de- 
pendent interactions. 

In this article, we investigate a quasi-2D Bose gas at 
zero temperature using the Bogoliubov theory of Bose- 
Einstein condensation. While the formalism is general 
to interactions of any long range or anisotropic form, we 
focus on the specific case where the dipoles are aligned 



perpendicular to the plane of kinetic freedom [125 . We 
note the formation of a roton in the excitation spectrum, 
due to a competition between the trapping potential and 
the attractive interaction in that direction. This roton 
introduces an instability in the system (the energy of the 
roton mode becomes equal to the energy of the ground 
state) at a critical value of the parameter; add^^^, where 
add is the dipole length scale (as defined in Ref. [8]), n 
is the projected 2D density (that is, the 3D particle dis- 
tribution projected down onto two dimensions), and Iz is 
the width of the density profile in the trapped direction. 
We pay particular interest to the implications this roton 
has on the quantum fluctuations in the system. We do so 
quantitatively by calculating the noncondensate fraction, 
the momentum distribution, the phase fluctuations, and 
the density fluctuations, as one approaches the critical 
point. Contrary to what might be expected from a the- 
ory which is essentially little more than an extension to 
mean field theory, we find intriguing predictions for the 
behaviour of the density fluctuations in the system near 
the critical point. The results suggest the emergence of 
a locally ordered state, in which density maxima lie on a 
hexagonal Br avals lattice. Although these predictions are 
qualitative, they give us clues regarding important con- 
cepts such as: How does the broken translational sym- 
metry in a solid emerge from a quantum fluid? Such 
questions regarding symmetry breaking and pattern for- 
mation are of interest in a wide range of disciplines, from 
cosmology to biology. 

A significant amount of recent theoretical focus has 
emerged for dipole-dipole interactions in systems of ei- 
ther bosons or fermions. In Ref. the Bose gas in its 
normal state is considered, and a semiclassical theory of 
the system above the transition temperature was used to 
develop a qualitative picture of the stability boundaries 
for various trapping geometries. Similar work was done 
in Refs. [m", Ts] for fermionic gases. At temperatures 
near quantum degeneracy, the question of Cooper-pair 
formation, and superfluidity, in dipolar fermionic sys- 
tems was considered in Ref. [l6|, [l7| . 



Inspired by tra- 
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ditional condensed matter achievements, it was theoreti- 
cally shown how one may obtain fractional quantum Hall 
states, and realisations of the Laughlin wavefunctions, in 
rapidly rotating dipolar Fermi gases [18]. Progress on 
low temperature (condensed) Bose systems has also been 
forth-coming, with a great deal of focus on numerical 
methods for efficiently solving nonlocal, nonlinear par- 
tial differential equations [19, 20]. Such equations ap- 
pear in the theory of the zero-temperature analysis of 
the condensate wavefunction [8]. Numerical solutions to 
the Bogoliubov-de Gennes equations in a trapped 3D gas 
were also studied Such solutions yield a direct un- 
derstanding of the excitation spectra and quantum de- 
pletion of such condensates. A point of interest in this 
work was the discovery of a roton-maxon in the excita- 
tion spectrum [22, 2^3] , as found in the neutron-scattering 
experiments with superfluid helium According to 

Landau's theory of superfluidity, this roton mode has a 
profound effect on the superfluid critical velocity of the 
system. This effect was studied in detail in Ref. ^ us- 
ing numerical solutions to the dipolar- Gross-Pitaevskii 
equation, with an object moving through the fluid. The 
results provide a powerful numerical validation of Lan- 
dau's theory. The anisotropy of this critical velocity was 
studied in Ref. [26|. Topological excitations of the con- 
densate mode have also been studied j27M3Q]. This yields 
the structure and stability of vortices and solitons in such 
sytems. 

The problem of stron gly correlated dipolar systems 
is a challenging one (3ll, l32|. It has long been under- 
stood, within the theory of classical 2D molecular dy- 
namics, that interactions of the form U{r) = D/r^ ex- 
hibit thermal melting transitions when the parameter 
F = D/{k^TrQ) (where ro is the average interparticle 
separation, T is temperature, and /cb is Boltzmann's con- 
stant) drops below a critical value, approximately given 
by Fc 62 ± 3 [33] . More recently, it was pointed out 
that, at low temperatures, quantum fluctuations can also 
melt the system [s^, [35|. When the ratio, r = Ei/Ek^ 
of average interaction energy; Ei = D/r^^ to zero-point 
kinetic energy; = /{mro) (where m is the mass of 
the particle), drops below the critical value; Tc — 18 ± 4. 
The results of both the molecular dynamics simulations 
at finite temperature, and the path-integral Monte Carlo 
method at zero temperature, indicate that the transition 
is first order, which suggests the possibility of a mixed 
phase regime. 

Both the classical and quantum phase transitions occur 
when the average interaction energy is large compared 
with the average kinetic energy. For this reason, theoret- 
ical attempts to understand the phase-transition rely on 
computationally intensive algorithms such as diffusion- 
Monte-Carlo, or path-integral-Monte-Carlo. While these 
calculations provide an exact treatment, we demonstrate 
that the considerably simpler Bogoliubov theory shows 
definite presursors to the transition, where a lattice struc- 
ture begins to emerge in the many-body density correla- 
tions, inspite of a uniform density remaining in the en- 



semble average. 

Our paper is organised as follows: In Sec. [TTlwe intro- 
duce the formalism and notation which we use through- 
out the manuscript. The formalism can be described 
as an extension of Bogoliubov's theory of Bose-Einstein 
condensation to account for long range and possibly 
anisotropic interactions. In Sec. lIIII we analyse the excita- 
tion spectrum of the system, paying particular attention 
to the formation of the roton. We identify the single pa- 
rameter which quantifies the strength of the roton in the 
system and the critical point at which the roton touches 
zero, and the system becomes unstable. In Sec. |IVl we 
study the noncondensate fraction as a function of the 
roton strength. We find this noncondensate fraction di- 
verges as one approaches the critical point. In Sec. |Vl 
we study the momentum distribution of the system, and 
find a rapid increase in the population of particles oc- 
cupying the roton modes as one approaches the critical 
point. In Sec. I VII we study the phase coherence as a func- 
tion of roton strength. In Sec. IVIII we study the density 
fluctuations as a function of roton strength. Finally, we 
conclude in Sec. IVIIII 



II. FORMALISM 

We start with the three-dimensional Hamiltonian for a 
bosonic field interacting via dipole-dipole interactions 
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where l^nt(r) = ^ / d^r''^\r')U{r - r')^(r') describes 
the two-body interactions. The field operator ^ sat- 
isfies the usual, equal-time, bosonic commutation rela- 
tions; [^(r),^t(r') = J(r - r'), and ^(r),^(r') = 

0. The two-body interaction potential 



^t(r),^t(r') 
is given by 
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where D is the dipolar coupling constant, and we have 
assumed that some external field aligns all dipoles along 
the z axis. The equation of motion for the Bose field op- 
erator is found from the Heisenberg equation of motion. 



(3) 



This equation is greatly simplified if we assume the 
field operator can be factorized into axial (along the z 
axis), and radial (in the x-y plane) components; ^(r) = 
?/^(p)(/)o(^), where p = {x^y). The form of ^o, in gen- 
eral will depend on the strength of the confining poten- 
tial ujz. For instance, if the confinement is very weak, 
such that huz <C /i (where /i is the chemical potential of 
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the system) then (/)o is weh approximated by an inverted 
parabola [sgI, Ist) of width Itf = \/2/i/ma;^ (the Thomas- 
Fermi profile). In the opposite limit, where huz ^ /i 
the system is purely 2D, and 0o is well approximated 
by the ground state harmonic oscillator wavefunction, 
with width /ho = yhjrmjjz. In this manuscript, for 
simplicity, we assume is given by a Gaussian pro- 
file, of width 1^. In general this value of 1^ can be 
determined by minimising the energy, using the ansatz 
^o(^) = (J^z\f^^~^^^ ^^^^7 and depending on the ra- 
tio hWz/lJi^ one will find /ho ^ Iz ^ Wf- We expect 
this theory to be quantitatively accurate in the 2D limit 
hwz ^ /i, yet still provide qualitatively correct predic- 
tions in the quasi2D and Thomas-Fermi regimes. 

Averaging Eq. (|3|) over the axial degree of freedom, we 
find 
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where V^'^^ip) = \ j p'^\p')U2j,{p - p')^{p'). The 
two body interactions are now described by the poten- 
tial U2'd{p) which includes an interaction with the mean 
density along the z-axis, and is given by 
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where h{x) = xe^ erfc(x). This interaction poten- 
tial asymptotes toward U2b{p) -> Dir^^'^ /{12V2) {Iz/pf 
when p ^ Iz^ but then becomes negative (attractive in- 
teractions) when p ^ Iz- This attractive component is 
directly due to the finite density of particles along the z 
axis. 

At zero temperature, we assume the single particle 
reduced density matrix has a macroscopically occupied 
eigenstate, which can be treated as a classical field (in 
an ensemble in which the phase-symmetry of ip has been 
broken). Under this assumption, the field operator can 
be split up into 



i,{p) = <t>{p)+Kp) 



(6) 



where (j) describes the macroscopically occupied state (the 
condensate wave function) and 5 describes all other states 
(which are assumed to have small occupation numbers). 
As a first order approximation, only the particles in the 
condensate are considered, and Eq. (j4]) reduces to the 
familiar dipolar Gross-Pitaevskii equation. 
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where uq = is the condensate density. 

The equilibrium solution to Eq. ^ \^ (j) ■ 
where the chemical potential is given by 
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and no denotes the projected 2D density (that is, the 
density of all particles projected onto the x-y plane). 

The noncondensate particles, which occupy states con- 
tained in 5 can be solved for, within the Bogoliubov ap- 



proximation, by substituting = 
Eq. dlj. Ignoring terms of order (5^ and higher (which 
describe interactions between noncondensate particles), 
and employing the Bogoliubov transformation 
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we obtain the following coupled equations for the quasi- 
particle amplitudes, 

EkUk = Eluk + no {uk + Vk) Tk [^2d] (10) 
-^kVk = ^ivk + ^0 {uk + Vk) Tk ^2-d\ (11) 

where E^, = -^^^ is the energy of a free particle and 
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, /— , 2 - 3y^h ( % ) is the Fourier trans- 
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form of the two-body interaction in Eq. (|5|). Equa- 
tions ([To]) and (pTj) can be solved, to give 
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and the excitation spectrum is given by 
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Using these solutions to the Bogoliubov-de Gennes equa- 
tions, we can investigate equilibrium correlation func- 
tions of the system, which is the aim of this paper. 



III. EXCITATION SPECTRUM AND 
STABILITY CRITERIA 

A particularly interesting feature of this system is 
the emergence of a roton in the excitation spectrum, 
Eq. (p!4|) . This spectrum can be rewritten by factoris- 
ing out the chemical potential and defining q = 
where ^ = h/ ^Jmp defines a healing length (or corre- 
lation length) for the fluid. This yields 
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where a = /^/^ is the ratio of confinement length to heal- 
ing length. So we find that the strength of the roton fea- 
ture in this system depends only on one parameter; a. 
Using the definition of ^ and p we can rewrite 
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FIG. 1: The excitation spectrum for quasi2D dipolar Bose 
gas given by Eq. ([T4j). 



FIG. 2: The noncondensate density for a quasi2D dipolar 
Bose gas as a function of a. The plot is generated from a 
numerical solution to Eq. (|18|) . 



where add = mD/ (121111?) is the dipolar length scale [11 . 
From Eq. ([16]) we see that the prominence of the roton in 
the spectrum depends on the number of particles within 
one dipole- length-scale ~ \/nadd^ multiplied by the num- 
ber of particles within one confinement width ~ yjnl^. 

We find that the short-range attractive component of 
[/2D renders the system unstable beyond a critical value 
of (J which we numerically determine to be <Jc ~ 1.8562. 
For G > Gc there exists complex values of E]^^ indicating 
a dynamical instability. 

Interestingly, the excitation spectrum at cr = is ex- 
actly equivalent to a system with contact interactions. 



IV. NONCONDENSATE DENSITY 

Calculating the noncondensate density provides an im- 
portant retrospective analysis of how accurate the as- 
sumptions of Bogoliubov theory are. The noncondensate 
density is given by, 
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and if we use q = ^k, and switch to polar coordinates, 
we find. 
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A necessary criteria to ensure our formalism is accu- 
rate, is that the noncondensate density be much less 
than the condensate density, n' /no <^ 1. Immediately 
from Eq. ([18]) we see that this criteria will be satisfied if 
^^no is sufficiently large. That is, the number of particles 
within one square healing length, needs to be large. Of 
course the ratio also depends on the value of the integral 
in Eq. (fT8|) . which depends on a. We numerically solved 



this integral and the result is shown in Fig.[2l We find the 
noncondensate density n' = l/47r when a = 0, and ini- 
tially decreases as a increases. However, when a becomes 
larger than ^ 1.5, the noncondensate density begins to 
rapidly diverge; ^ oo as a ^ dc. This divergence is 
due to a rapid increase in the population of the roton 
modes. This fact will become clear after examining the 
momentum distribution of the system. 



V. MOMENTUM DISTRIBUTION 

We calculate the momentum distribution of the system 
as a means of analysing the occupation number of the 
roton modes. In general, this can be found via the Fourier 
transform of the one-body density operator. 
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where A is the total area of the 2D system. Within our 
current formalism, this reduces down to 
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The first term in Eq. (|2Q|) is simply the condensate mode, 
and the second term constitutes the so-called quantum 
depletion caused by the interactions. This quantum de- 
pletion is zero, only in the case where there are no inter- 
actions in the system and Ek = E^. 

We analyse the asymptotic behaviour of 77 (k), and find 
that for /c^ <C 1, the distribution depends weakly on the 
value of cr; 



r]{k) ^ no6{k) 
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However, in the region 0.75 < /c^ < 1 where the roton ap- 
pears in the spectrum, the momentum distribution shows 
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FIG. 3: The momentum distribution, Eq. (|2Qp , for a quasi2D 
dipolar Bose gas for different values of a. The distribution 
shows a remarkable sensitivity to a as it approaches Gc- 



FIG. 4: The first order correlation function, g'^^^ of Eq. ([21]) . 
for a quasi2D dipolar Bose gas for different values of a. The 
phase coherence shows an unusual oscillatory behaviour as 
(J ac. 



a remarkable sensitivity to the value of a, as shown in 
Fig. m As cr approaches ctc, we see a dramatic increase 
in the population of these roton modes. This behaviour 
is similar to that shown in a quasi- ID geometry for laser 
induced dipoles [38] . This increase in population explains 
the rapid divergence of the noncondensate density seen in 
Fig. [21 As the system approaches the instability, macro- 
scopic occupations of these roton modes begin to exist, 
resulting in the systems eventual collapse. 



VI. PHASE COHERENCE 



In this section, we calculate the first order correlation 
function for the system 



(22) 



where n = no + nMs the total density of the system. 
This quantity tells us how correlated the phase of a 
given particles wavefunction at some point p is with the 
phase of that same single-particle-wavefunction at an- 
other point p^ A Bose-Einstein condensate is said to 
exist when g^^\p) decays to some finite (nonzero) value 
as p ^ oo. A pure Bose-Einstein condensate would imply 
g^^\p) = 1 everywhere, exactly as in Glauber's original 
definition of a coherent state. Interactions in the system 
destroy the coherence, in our case this is given by 
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Again we switch to the dimensionless variable q 
and convert to polar coordinates, to find 

1 



no 
n 



2 ' 



2q 



(24) 

where Jo is a Bessel function. For the case a = 0, where 
the spectrum resembles that of contact interactions, the 
phase coherence is given by, 

1 
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where Ii and Ki are Bessel functions. For a 7^ 0, we 
solve the integral in Eq. (|24|) numerically, and plot the 
results in Fig. HI Interestingly, the correlation function 
displays an oscillatory nature for values of a close to ac. 
This oscillation implies partial-collapse/revival type be- 
haviour, as if the system is attempting to fragment into 
smaller pieces. 



VII. DENSITY FLUCTUATIONS 

A. Two-body correlations 

Two-body, density-density correlations yield impor- 
tant information with respect to the linear-response of 
the system to external perturbations. The quantity of 
interest to us, is the so-called g^'^\p^ p') function, which 
quantifies the conditional probability of detecting a par- 
ticle at p, given that one has been detected at p' . This 
quantity is defined as 
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which in our case reduces down to 



(27) 

Switching to the dimensionless variable q = k^, and con- 
verting to polar coordinates, we find 
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The integral in Eq. (|28]) converges for all p 7^ 0. When 
p = however, the integral has an ultraviolet diver- 
gence. The physical reason for this divergence is well 
known, and is related to the incorrect expression for the 
short-range interaction in U^d- Physically, the interac- 
tion potential [/2D should include a short-range van der 
Waals type interaction, which would introduce a short- 
wavelength cut-off (a finite upper-Hmit) to the integral 
of Eq. ([28|) . This ultraviolet cut-off does not significantly 
alter the behaviour for p ^ 0, and the local correlation at 
p = is only logarithmically dependent on this cut-off. 

The numerical solutions to the integral in Eq. (|28]) are 
shown in Fig. O We see that, when a = 0, the case 
is equivalent to repulsive contact interactions, where the 
repulsion between particles causes a decrease in the prob- 
ability of detecting two particles close to one another i.e. 
g^'^^ < 1. As cr increases, particles begin to feel the ef- 
fect of the third dimension, that is, the attractive core 
in the interparticle potential /72d- This attraction causes 
an increase in the local value of g^'^^ > 1. The most in- 
triguing behaviour in g'^'^^ appears near the critical point, 
a ^ (Tc. Near this point we see g^'^\p) dip below unity, 
as one would expect for repulsive long range interactions, 
but then oscillate between local maxima/minima which 
are above and below unity respectively. The first and 
highest peak occurs at p = 0, and is caused by the rem- 
nant attractive component of the 3D interactions. The 
second peak, occurs at a distance rmax ^ 7.92^, and can 
be loosely construed as a preferred separation distance 
between particles (or possibly clumps of particles). Sub- 
sequent peaks occur at approximately integer multiples 
of Tmax- This bchaviour seems qualitatively similar to 
what one might expect from a fluid which is undergoing 
a transition to a solid state, where long range diagonal 
order would be present. Because the translational sym- 
metry of the system remains unbroken, the diagonal or- 
der becomes off-diagonal order in the g^'^\p^ p') function. 
The order which we find in this case, however, is not long 
range, and decays proportional to l/\p — p' 




FIG. 5: The second order, density-density correlation func- 
tion, g^'^^ of Eq. ()28p . for a quasi-2D dipolar Bose gas for 
different values of a. The different colors correspond to dif- 
ferent values of a, in (a) the black line is a = 0, red line is 
a — 0.38, green line is a = 0.75, and blue line is a = 1.5. In 
(b), the black line is cr = 1.81, and red line is a = 1.8562. 



B. Higher order correlations 

To further test the hypothesis that the system is at- 
tempting to form an ordered state, we also study higher 
order correlations in the system, three-body, four-body, 
etc. In general the Nih order density correlation is given 

by, 



9^ '{Pi,---,Pn) = , 

(29) 

this can be interpreted as a conditional probability of 
detecting a particle at position p^, given the detection 
of N — 1 particles at positions • • • 7 Pat- Within our 
current theory, this function can be simplified down to 
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where the notation; j) — ^^=2 Xli=i summa- 
tion over all pairs of indices, and the arguments of g^^^ 
have been omitted for notational convenience. 

To visualise this high dimensional observable, we 
choose fixed values for p^, up to p^v-i follow- 
ing way. The first particle is detected at some arbitrary 
point. This choice is arbitrary, as there is complete trans- 
lational symmetry in the system at this point. Subse- 
quent particles are then chosen at the most likely points 
given by the behaviour of g^'^\pi,P2), g^^\p^, P2, p^), 
etc, with the proviso that these new points are spatially 
separated from the previous points. For instance, con- 
sider Fig[6ja), where we plot g^^\pi^ p2^ P)- We begin by 
(arbitrarily) setting p^ = (— rmax/2,0). We know then, 
from our calculation of g^'^'^ in the previous subsection, 
that there will be an increased probability of detecting 
a second particle somewhere on the ring centred at pi, 
and with a radius of rmax [see the red line in Fig. Efb)]. 
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So we choose to place the second particle at the point 
p2 — (^max/2,0), such that the two particles are sepa- 
rated by a distance of exactly rmax- We then plot the 
function g^^\pij p2-, p) — 9^'^\pi^ P2) ^^r this choice of p^ 
and shown in Fig.[6fa) (note that we have subtracted 
off the lower-order g'''^\pn P2) so that the asymptotic 
value of the plotted function ^ 1 as p ^ 00). The black 
crosses in the figure indicate the positions of p^ and p2- 
As was the case, with g'^'^\ we again see a series of max- 
ima/minima which are above and below zero respectively. 
The maxima along the x-axis are placed at integer mul- 
tiples of (rmax/2,0). However, along the ?/-axis the first 
maxima is at approximately (0, ir^ax a/3/2), consistent 
with the formation of a triangular lattice (the densest 
packing arrangement of spheres in 2D). 

Continuing on in this fashion, we construct higher or- 
der correlation functions, and visualise them in this way. 
The results are shown in Fig. [6] for a variety of different 
choices of p^^, . . . , Ptv-i- The results suggest the forma- 
tion of localised order consistent with a triangular lattice 
of high-density droplets, separated by a distance of ap- 
proximately Tmax- The value of a is very close to the crit- 
ical value (Jc, where the theory becomes invalid. For this 
reason, we proceed with a degree of scepticism, and only 
claim that, although the quantitative aspect of these pre- 
dictions is questionable, the qualitative aspect seems very 
reasonable. Specifically, the results give us clues as to the 
ground state of the system past the phonon-instability. 
This problem was studied numerically in Ref. [39^ , where 
they conluded that, beyond the critical point, the sys- 
tem forms a soliton gas. These conclusions were reached 
by simulating the dipolar Gross-Pitaevskii equation, with 
an instantaneous quench in the dipolar coupling constant 
such that a > ac. We conject that it may be possible to, 
instead, form a state in which the solitons spatial distri- 
bution is ordered into a triangular lattice, somewhat like 
a soliton solid. This state could be reached if the dipolar 
coupling constant is increased adiabatically, or at least 
very slowly, beyond <Jc. 

VIII. CONCLUSIONS AND DISCUSSION 

In this article, we have formulated a simple theoretical 
model for predicting the behaviour of quasi-2D dipolar 
Bosons close to zero temperature. We solved the model 
in the case of a uniform system in equilibrium, and stud- 
ied the excitation spectrum, the momentum distribution, 
the phase fluctuations and the density fluctuations. Re- 
garding the excitation spectrum, roton modes form in 
the vicinity of k ^ 0.87/^ due to a competition between 
the trapping potential and the attractive interaction in 
that direction. We found a single parameter, a defined in 
Eq. (p!6|) . which determines the strength of the roton, and 
the point of collapse. For the momentum distribution, we 
found it decays as 1/k when /c^ <C 1. However, as the 
system approaches the critical point, the occupation of 
the roton modes begin to increase rapidly. Eventually 



these roton modes become macroscopically occupied re- 
sulting the eventual collapse of the system. The level of 
phase fluctuations which occur in the system is compa- 
rable to the fraction of noncondensate particles in the 
system. This initially decreases with increasing a, but 
then rapidly diverges as the roton modes become highly 
occupied. The density fluctuations display a variety of 
different behaviour due to competition between the inter- 
atomic repulsion at large distances versus the short range 
attraction at short distances (due to the mean field of 
particles in the trapped direction). The interplay of these 
different length scales manifests as the behaviour shown 
in Fig.[5fa). As a approaches the critical point, the den- 
sity fluctuations of ^^^^ begin to show an anomalous os- 
cillatory behaviour, with maxima and minima above and 
below unity, respectively. We construe this oscillatory be- 
haviour as a preferred separation distance between parti- 
cles or clumps of particles. We further pursued this idea 
by calculating higher order density fluctuations, shown 
in Fig. [6l The coexistence of crystalline order and su- 
perfluidity (supersolidity) has been a subject of intense 
recent debate [40| . Although what we observe here is not 
true crystalline order, since it is not long range, it does 
provide an intriguing possibility for the manifestation of 
a supersolid state. 

A. Relation to experiments 

We conclude with some brief remarks on the possibil- 
ity of experimental observation of our results. To date, 
dipolar Bose-Einstein condensates have been achieved 
with Chromium [5], Dysprosium [6j, and Erbium 0, 
which have dipole length scales of approximately 0.85nm, 
7.7nm, and 3.9nm respectively. We assume 2D densities 
of around lO^^m"^ to lO^^m"^. These estimates are ap- 
proximately based on what has been achievable in tradi- 
tional cold atom systems such as Rubidium We see 
then that the if the confinement width, is on the order 
of 2-3/im then a will lie in the stable regime, and by tun- 
ing Iz one should be able to observe the roton facilitated 
collapse. 

A viable alternative to ultra-cold gases is in systems of 
strongly correlated electrons, such as the exciton gases 
found in certain semiconducting heterostructures [42|. 
Such systems have been studied using a path-integral 
Monte-Carlo method in Ref. (isf . 
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FIG. 6: Higher order density fluctuations for a — 1.8562. In (a)-(j) we plot g^^\pi, . . . ,pj^_i,p) — g^^~^\p'^^ . . . ,Pn_i), for a 
variety of different values of N, aswell as different choices of p^^, . . . , pjv-i- The black crosses show the positions of p^^, . . . , p^v-i? 
which are chosen according to the procedure described in the text. From the number of black crosses one can deduce the value 
of N in each subfigure. We note the intriguing emergence of hexagonal order in the density distributions. 




FIG. 7: Higher order density fluctuations for a = 1.8562, as 
in Fig. (S^e) shown as a surface plot with contour lines. 
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